Invasion percolation: new algorithms and universality classes

Employing highly efficient algorithms for simulating invasion percolation (IP), whose execution time scales as O[M log(M)] or better for a cluster of M sites, and for determining the backbone of the cluster, we obtain precise estimates for the fractal dimensions of the samplespanning cluster, the backbone, and the minimal path in order to identify the universality classes of four different IP processes (site and bond IP, with and without trapping). In two dimensions IP is characterized by two universality classes, one each for IP without trapping, and site and bond IP with trapping. In a three-dimensional site IP with and without trapping is in the universality class of random percolation, while bond IP with trapping is in a distinct universality class, which may be the same as that of optimal paths in strongly disordered media. Multiphase flow phenomena in porous media are relevant to many problems of great scientific and industrial importance, including extraction of oil, gas and geothermal energy from underground reservoirs, food and soil sciences, powder technology and materials science [1]. Invasion percolation (IP), a model introduced [2] for describing the evolution of the interface between an invading and a defending fluid in a porous medium, has provided deep insight into such phenomena. In addition, IP is relevant to a host of other problems, including characterization of optimal paths and domain walls in strongly disordered media [3, 4], and even simulation of the Ising model at the critical temperature [5]. Moreover, IP is one of the simplest parameter-free models which exhibits self-organized criticality [6], another subject of current interest. Two main variants of IP have been studied so far. In nontrapping IP (NTIP) the defending fluid is compressible and the invading fluid can, potentially, enter any region on the interface which is occupied by the defending fluid. In the second and more common case, the defending fluid is incompressible and is trapped if a portion of it is surrounded by the invading fluid. In addition to the compressibility, one must also take into account the capacity of the fluids to wet the internal surface of the medium [1]. The process by which a wetting fluid is drawn spontaneously into a porous medium is called imbibition, while the forcing of a nonwetting fluid into the pore space is called drainage. We model the porous medium as a network of pores or sites connected by throats or bonds which have smaller radii than the pores. In IP, the potential displacement events are ranked according to the capillary pressure threshold that must be exceeded before that event takes place. During imbibition, the invading fluid is drawn 0305-4470/99/490521+09$30.00 © 1999 IOP Publishing Ltd L521 L522 Letter to the Editor first into the smallest constrictions, for which the capillary pressure is large and negative, and it goes last into the widest pores. Displacement events are therefore ranked in terms of the largest opening that the invading fluid must travel through, since it is from these larger capillaries that it is most difficult to displace the defender. Imbibition is therefore a site IP and, in contrast, drainage in which the invader has most difficulty with the smallest constrictions is a bond IP. Important differences arise in the structure of the invading fluid paths depending on whether one considers NTIP or TIP [2, 7, 8]. The question of the universality class of IP has not conclusively been established, and the literature contains contradictory claims about it. The scaling properties of NTIP are believed to be consistent with random percolation (RP). For trapping IP (TIP) the fractal dimension Df of the sample-spanning cluster (SSC) in two dimensions is smaller than that of RP [2]. In three dimensions no significant difference between TIP and RP has been reported. It was originally argued [2] that the fractal properties of IP do not depend on whether one simulates a site or bond IP. Recently, however, it was argued [9] that important differences arise in the structure of the invading fluid’s paths when comparing site and bond IP. Porto et al [4] used a mapping from the minimal (shortest) paths in TIP to the optimal paths in strongly disordered media, and presented numerical evidence that for TIP the fractal dimensionDmin of the minimal path in all dimensions is not the same as that of RP, and hence they argued that TIP and RP do not belong to the same universality class. Barabási [10] argued that loopless bond TIP is in the universality class of RP. The reason that the universality class of IP has not conclusively been established is the relatively limited accuracy of the numerical results reported so far. This is due to the difficulty of simulating IP models, especially TIP, with very large lattices, since for all the present IP algorithms the simulation time grows as O(NM), where N is the number of sites (or bonds) in the lattice and M the number of sites in the cluster. This has severely hampered studies of IP models. Since the differences between values of various fractal dimensions and scaling exponents of IP and RP appear to be small, it is critical to be able to use very large lattices to establish the universality classes of IP models. The purpose of this letter is twofold. First, we describe a highly efficient algorithm for simulating IP in which the simulation time grows as O[M log(M)], which enables us to use very large lattices. Second, using the new algorithm we measure the fractal properties of various IP models with an accuracy which is at least an order of magnitude better than the best previous estimates. This enables us to test various hypotheses and conclusively establish the true universality classes of the IP models. In the conventional algorithms [1, 2, 4, 8] the search for the trapped regions is done after every invasion event using a Hoshen–Kopelman algorithm [7, 8], which traverses the whole lattice, labels all the connected regions, and then only those sites (bonds) that are connected to the outlet face are considered as potential invasion sites (bonds). A second sweep of the lattice is then done to determine which of the potential sites is to be invaded in the next time step. Thus, each invasion event demands O(N2) calculations. This is highly inefficient for two reasons. First, after each invasion event only a small local change is made to the interface; implementing the global Hoshen–Kopelman search is unnecessary. Secondly, it is wasteful to traverse the entire network at each step to find the most favourable site on the interface since the interface is largely static. We tackle the first problem by searching the neighbours of each newly invaded site to check for trapping. This is ruled out in almost all instances. If trapping is possible, then several simultaneous breadth-first ‘forest-fire’ searches are used to update the cluster labelling as necessary [11]. This restricts the changes to the most local region possible. Since each site can be invaded or trapped at most once during an invasion, this part of the algorithm scales as O(N). This cluster searching method has some similarities with the ‘perimeter scouting’ algorithm for two-dimensional (2D) clusters [12, 13]. In this algorithm one checks whether the most recently invaded sites could have been trapped in the interior of Letter to the Editor L523 the cluster. If so, oriented walks are started on the just invaded site, pointing away from it to the neighbouring sites which are those that neither belong to the cluster nor are candidates for invasion. The walks continue until all but one of them have again reached the site of origin. The growth sites visited by these walks are then eliminated from the list of active sites. This method is effective in two dimensions but not as efficient in three dimensions. Our method differs from it by searching cluster volumes rather than perimeters, and incorporating local checking to minimize cluster searching and is thus equally effective in three dimensions. The second problem is solved by storing the sites (bonds) on the fluid–fluid interface in a list, sorted according to the capillary pressure threshold needed to invade them. This list is implemented via a balanced binary search tree, so that insertion and deletion operations on the list can be performed in log(n) time, where n is the list size. Sites (bonds) that are designated trapped using the above-described procedures are removed from the invasion list. Each site is added and removed from the interface list at most once, limiting the cost of this part of the algorithm to O[N log(n)]. Thus, the execution time forN sites (bonds) is dominated (for large N ) by list manipulation and scales at most as O[N log(N)]. While the execution time is approximately O[N log(n)], in practice the time and memory requirements depend on the total number of lattice sites and those forming the cluster boundary. For example, we find empirically that for 3D TIP the execution time scales as M1.24, and the memory use is 20 bytes for each lattice site plus 64 bytes for each cluster site. On a 500 MHz 21164A Alpha microprocessor, a trapping cluster of 2×105 sites is grown on a 181×181×181 lattice in 12.0 s, using 120 Mb of memory, while in two dimensions a cluster of 5×105 sites is grown on a 2000×2000 lattice in 12.0 s, using only 52 Mb. Complete details of the algorithm, which can be used for arbitrary networks, are given elsewhere [14]. We have also used a new optimized algorithm to identify the minimal path length, the sites comprising both the elastic backbone [15], i.e., the set of the sites that lie on the union of all the shortest paths between two widely separated points, and the usual transport backbone, so that the backbone search and computations do not affect the overall execution time of the algorithm. In the past, numerous algorithms have been proposed in the literature [15–18], most of which are too slow or limited to 2D systems. For example, a recent method [19] that uses a matching algorithm takes longer to identify the backbone than the percolation algorithm used here takes to generate it. An alternative method was recently presented by Babalievski [11], based on depth-first searching out from the elastic backbone [15] to identify loops of occupied sites. This method works well for low connectivity clusters but loses efficiency where the SSC is composed of large well-connected regions. In the latter case, some sites need to be visited numerous times before their status is decided. The method used here is an optimization of this in which the distance on the cluster from the inlet face to each cluster site is used to guide the depth-first search. In this algorithm, there are three major steps which are as follows. (i) Using a breadth-first search algorithm, we label each site in the cluster with its ‘cluster distance’ from the inlet face, and then use this information to burn backwards from the outlet face and identify the elastic backbone. At the same time, we construct the ‘branch points list’—a list of all the cluster sites that are adjacent to the elastic backbone but are not part of it. The branch points list should be ordered with the sites closest to the inlet face listed first. Note that the elastic backbone sites are part of the backbone. (ii) We stop if the branch points list is empty. Otherwise, we perform a depth-first search from the last site in the branch points list, flagging all the sites that are visited. During the search, unexplored branch points are added to the branch points list, while another list tracks the sites that have been flagged as visited. We then perform an important optimization during the depth-first search: if there are multiple branches from a single site, the site labelled as L524 Letter to the Editor being closest to the inlet face is always the first to be explored. (iii) The depth-first search terminates when one of two conditions are satisfied: (1) the search contacts the backbone again at a different site from whence it started, in which case the sites in the visited-sites list are flagged as backbone sites; or (2) it retreats back to its starting site, at which point there will be no sites left in the visited-sites list. (iv) The algorithm continues at step (ii). In this way the elastic backbone, the transport backbone and the dangling ends of the SSC are all identified. Examples of execution times for this algorithm running in three dimensions on a 533 MHz 21164A Alpha processor are 0.02 and 0.12 CPU s for 323 and 643 lattices, respectively. The cluster on which these calculations were performed was a SSC generated by a NTIP. When compared with the timings reported in [19] on an equivalent hardware, our algorithm is faster by a factor of 7 for the 323 lattice and by a factor of 12 for a 1283 lattice, and thus the larger the lattice size, the more efficient is this algorithm. Comparison with the other algorithms [15–17] is also favourable. We used Ld−1 × 2L lattices in d dimensions with reflecting boundary conditions on the edges. Cluster properties were measured within the central L region. We also considered lattices of size Ld−1 × nL in d dimensions, where n = 4, 8, 16, and measured the cluster properties within the central L region to check for end effects. We also studied the effect of periodic versus reflecting boundary conditions. The effect on the results was negligible. The simulations reported in this letter used up to 16 3842 and 5123 lattices in two and three dimensions, respectively, larger than the largest previously-used lattices by a factor of 5–10. At least 2000 realizations were used for the largest lattices and over 105 realizations for smaller L. We first consider the fractal dimension Df of the SSC. If we define a local fractal dimensions Df (M) ≡ d lnM/d lnL, then according to finite-size scaling (FSS) Df (M) will converge to the asymptotic value for large M according to, |Df −Df (M)| ∼ M−α , where α is a correction-to-scaling exponent to be estimated from the data. Combining the definition of Df (M)with the FSS equation gives us a differential equation which has an analytical solution: c1 +DfM α = c2L (1) where c1 and c2 are constants. We then fit the data to equation (1) to estimate both Df and α simultaneously. By so doing we also avoid statistical pitfalls of the two-stage process used in [13] in which the data are first divided into various bins, the local Df (M) are estimated by numerical differentiation of lnM with respect to lnL, and then Df (M) is plotted versus M−α and α is varied until the the ‘best’ straight line is obtained. In addition, we also obtain reliable estimates for the confidence intervals of α and the fractal dimensions (see below). From this analysis we obtained, Df ' 1.8949 ± 0.0009 (α ' 0.475), completely consistent with the exact value,Df = 91 48 = 1.895 83 for RP. For TIP we obtained,Df ' 1.825± 0.004, more than an order of magnitude more precise than the previous estimates [2–4]. In three dimensions we found that, as far asDf is concerned, there is no significant difference between NTIP and TIP, hence confirming explicitly that trapping occurs rarely in three dimensions; we found Df ' 2.528 ± 0.002, consistent with the most recent estimate of Df for RP [20], Df ' 2.523± 0.004. Moreover, no difference was found between values of Df for site and bond TIP. However, the universality classes of IP are not determined solely by Df . The fractal properties of the backbone and the minimal paths, which are particularly important for flow and transport processes, must also be considered. Our simulations indicate that, as a consequence of the trapping and the choice of site or bond IP, strong differences exist between the backbone Letter to the Editor L525 Table 1. Fractal dimensions for invasion percolation in two dimensions. Numbers in parenthesis are the estimates obtained from local Df (M) analysis.